Moving body observation method

ABSTRACT

A moving body observation method includes acquiring time-series images taken at a predetermined time interval, obtaining an estimated position of a moving body in the time-series image based on an estimated motion of the moving body, extracting an image including the estimated position of the moving body from images of the time-series images, as an extracted image, generating a template image by stacking the extracted images in such a manner that reference points of the moving body in the extracted images coincide with each other, and identifying a position of the moving body in each of the time-series images by performing template matching using the template image for at least one of the time-series images.

CROSS REFERENCE TO RELATED APPLICATIONS

This application is a continuation application of International Application No. PCT/JP2017/045949, now WO2018/230016, filed on Dec. 21, 2017, which claims priority to Japanese Patent Application No. 2017-115833, filed on Jun. 13, 2017, the entire contents of which are incorporated by reference herein.

BACKGROUND 1. Technical Field

The present disclosure relates to an observation method for a moving body affected by fluctuation or the like.

2. Description of the Related Art

Conventionally, when an image of an object (hereinafter, “moving body”) that moves under a predetermined physical law (for example, at a constant velocity) is taken as an observation object by a device such as a CCD camera, it is difficult to find the moving body because the moving body is buried in background noise, for example, if the moving body is at a far distance or is small, or if brightness of the moving body is weak. Particularly, because the velocity of the moving body is unknown, it is impossible to change the orientation of the camera in accordance with the moving body even if long exposure is performed. Therefore, a motion blur of the moving body occurs in a field of view, so that the image of the moving body cannot be taken.

Against the above problem, Japanese Patent Application Laid-open No. 2014-51211 (Patent Literature 1) proposes a moving body detection method that removes noise by using images acquired in time series. This method uses a number of images taken at a constant interval and performs the following processes for every sequence of consecutive images. (a) Assuming that a moving body performs a uniform velocity motion at a certain velocity, images are shifted at an assumed velocity and are stacked on each other, a slightly bright portion that is common to data at the same position is found, and an S/N ratio (a noise-to-signal ratio) is improved. (b) The assumed velocity is changed over the entire expected velocity range of the moving body, and the process (a) described above is repeated. A velocity and a position at which a bright point is found is acquired as a detection result. (c) From the results obtained by the processes (a) and (b) for every image sequence, those that are consistent with each other (for which the positions and the velocities are not discontinuous) are integrated as the same moving body.

SUMMARY

The moving body detection method in Patent Literature 1 described above uses a so-called stacking method, and can detect the position of a moving body by using a number of images. However, this method cannot recognize the locus or the behavior of the moving body. For example, even if N images are used in the stacking method, only the position of the moving body in an n-th (n<N) image can be detected, and the position of the moving body in each of the remaining (N−1) images cannot be detected. Therefore, use of the stacking method only provides very sparse information in terms of time.

The present disclosure has been made in view of the problems described above, and it is an object of the present disclosure to provide a moving body observation method that can accurately detect a moving body by using a stacking method.

An aspect of the present disclosure is a moving body observation method comprising: acquiring time-series images taken at a predetermined time interval; obtaining an estimated position of a moving body in the time-series images based on an estimated motion of the moving body; extracting images each including the estimated position of the moving body from a plurality of images of the time-series images, as an extracted image; generating a template image by stacking the extracted images in such a manner that reference points of the moving body in the extracted images coincide with each other; and identifying a position of the moving body in at least one image of the time-series images performing template matching using the template image for the at least one image of the time-series images.

Each of the reference points of the moving body in stacking of the extracted images may be the estimated position of the moving body.

Each pixel value in the template image may be set to any one of an average pixel value, a median, and a mode of a corresponding pixel in the extracted images.

The moving body observation method may further comprise binarizing each pixel value in the time-series images or the extracted images; and calculating a center of gravity of the moving body in each of the binarized time-series images or each of the binarized extracted images. In this case, each of the reference points of the moving body in stacking of the extracted images may be the center of gravity of the moving body.

The moving body observation method may further comprise binarizing each pixel value in the time-series images or the extracted images; and calculating a maximum-brightness position of the moving body by fitting a normal function to two one-dimensional distributions obtained by integrating a pixel distribution of the moving body in each of the binarized extracted images in predetermined two directions. In this case, each of the reference points of the moving body in stacking of the extracted images may be the maximum-brightness position of the moving body.

The moving body may be debris on an earth orbit or a microorganism in a liquid.

According to the present disclosure, it is possible to provide a moving body observation method that can accurately detect a moving body by using a stacking method.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a block diagram of a moving body observation device according to an embodiment of the present disclosure.

FIG. 2 is a flowchart of a moving body observation method according to the embodiment.

FIG. 3 is an explanatory diagram of each step of the moving body observation method according to the embodiment.

FIGS. 4A and 4B are respectively an example of an intensity distribution in an image obtained by the moving body observation method according to the embodiment, where FIG. 4A illustrates an intensity distribution in an extracted image, and FIG. 4B illustrates an intensity distribution in a template image Q_(T).

FIG. 5 is a flowchart according to a first modification of the embodiment.

FIG. 6 is a flowchart according to a second modification of the embodiment.

DESCRIPTION OF THE EMBODIMENTS

An embodiment of the present disclosure will be described below with reference to the accompanying drawings. FIG. 1 is a block diagram of a moving body observation device according to the present embodiment. As illustrated in FIG. 1, a moving body observation device 10 includes an observing device 20 and a controller 30. The observing device 20 includes an optical system 21 and an imaging unit 22. The optical system 21 forms an optical image including a moving body P that is an observation object, and is a telescope or a microscope, for example. The imaging unit 22 acquires the optical image generated by the optical system 21 as image data, and is a CCD camera, for example. The controller 30 includes a CPU (a processor) 31, a storage unit 32, and an input/output unit (I/O) 33. The CPU 31 performs various types of calculation using a program and data stored in the storage unit 32, for example. The storage unit 32 stores therein a program according to the present embodiment and various types of data, and is configured by a semiconductor memory, a hard disk drive (HDD), or the like. For example, the storage unit 32 also stores therein image data sent from the observing device 20. The input/output unit 33 performs communication of data, an instruction, and the like with an external device, and is connected to the imaging unit 22 of the observing device 20 in the present embodiment. Therefore, image data acquired by the imaging unit 22 is stored in the storage unit 32 via communication between the imaging unit 22 and the input/output unit 33 and is used for calculation by the CPU 31.

FIG. 2 is a flowchart of a moving body observation method according to the present embodiment. FIG. 3 is an explanatory diagram of each step of the moving body observation method according to the present embodiment. As illustrated in FIGS. 2 and 3, the moving body observation method according to the present embodiment includes an image acquiring step (Step S10), a position estimating step (Step S20), an image extracting step (Step S30), a template generating step (Step S40), and a position identifying step (Step S50).

The moving body P that is an observation object is debris on the earth orbit, for example. On the earth orbit, there are orbiting artificial satellites with various purposes, for example, a military satellite, a communications satellite, a science satellite, an observation satellite, and a navigation satellite. When these artificial satellites become no longer functional because of a failure, or finish their roles and end their service life, they are left on the orbit as they are to become debris (also referred to as “space dust” or “space debris”) that move on the orbit in many cases. Further, remains of rockets or the like used for launching artificial satellites or the like is also left on the orbit as debris.

Currently, the number of debris that is orbiting the earth reaches several thousands or more, and enters a self-multiplying stage in which the number of debris is increased by natural collision. It is likely that this debris collides with an artificial satellite that is being used when the satellite is orbiting or an artificial satellite that is being launched, and therefore an observation method with high accuracy is required.

The observation object in the present embodiment is not limited to the debris described above. That is, the moving body observation method can be applied to an object with fluctuation of a shape in images acquired at a predetermined time interval as the observation object. Therefore, the moving body P may be a microorganism in a liquid such as a chemical agent (for example, a still microorganism). In this case, a microscope is used as the optical system 21. The fluctuation is caused because a moving body is at a far distance, the size of the moving body is small, or the brightness of the moving body is weak, for example.

The image acquiring step (Step S10) acquires multiple (for example, N) time-series images F_(n) (n=1 to N) including the moving body P, which are taken at a predetermined time interval. In the present embodiment, the observing device 20 of the moving body observation device 10 takes multiple (for example, N) images including the moving body P at a predetermined time interval, and inputs image data of those images to the controller 30 (the storage unit 32) of the moving body observation device 10 as time-series images F₁ to F_(N).

After each time-series image F_(n) is acquired, a process of removing a bright pixel, such as a star, which may cause false detection, may be performed for each time-series image F_(n). For example, it is possible to remove a noise component of a star included in an image signal by translating the images in accordance with a moving direction and a moving amount of the star, stacking the images on each other, and subtracting pixel values of the same pixels that are stacked.

The position estimating step (Step S20) obtains estimated positions L_(n) (n=1 to N) of the moving body P in the time-series images F_(n) (n=1 to N) based on an estimated motion of the moving body P. For example, assuming that the motion of the moving body P is a uniform linear motion, the estimated motion can be represented by the following formula 1.

$\begin{matrix} \left\{ \begin{matrix} {{{PosX}(k)} = {{{PosX}(n)} + {\left( {k - n} \right) \times {VelX}}}} \\ {{{PosY}(k)} = {{{PosY}(n)} + {\left( {k - n} \right) \times {VelY}}}} \end{matrix} \right. & \left( {{Formula}\mspace{14mu} 1} \right) \end{matrix}$

Here, PosX(k) represents an X coordinate of an estimated position L_(k) of the moving body P in a time-series image F_(k) (a frame number k), PosY(k) represents a Y coordinate of the estimated position L_(k) of the moving body P in the time-series image F_(k), VelX represents an X-direction component of a velocity of the moving body P in the time-series image F_(n), VelY represents a Y-direction component of the velocity of the moving body P in the time-series image F_(n), and k represents a positive integer that is equal to or less than N and is not n. PosX(n) and PosY(n) are an X coordinate and a Y coordinate of the estimated position L_(n) of the moving body P in the time-series image F_(n), respectively, and are selected from any one of the acquired time-series images. The estimated position L_(n) may be referred to from an observation history of the moving body P stored in the storage unit 32 or may be calculated again together with an estimated moving velocity of the moving body P. In the latter case, for example, while the estimated moving velocity of the moving body P is sequentially changed in a range of (−VelX, −VelY) to (VelX, VelY), the images are stacked on each other. A velocity at a time when the degree of coincidence between candidates of the moving body P in the images is high is set as the estimated moving velocity, and the candidates are set as the moving body P.

Assuming that the moving body P performs a circular motion, its estimated motion is represented by the following formula 2 or 3. These circular motions each represent a circular orbit with a small radius of curvature that can be approximated by a line in either of the X direction and the Y direction. Here, A represents an X coordinate of the center of the circular orbit described above, and B represents a Y coordinate of the center of the circular orbit described above. Formula 2 represents a case where VelX≤VelY and Formula 3 represents a case where VelX>VelY.

$\begin{matrix} \left\{ \begin{matrix} {{{PosX}(k)} = {{{PosX}(n)} + {\left( {k - n} \right) \times {VelX}}}} \\ {{{PosY}(k)} = {{{PosY}(n)} + {\left( {k - n} \right) \times {VelX} \times}}} \\ {\left( {{{PosX}\left( {k - 1} \right)} - A} \right)/\left( {{{PosY}\left( {k - 1} \right)} - B} \right)} \end{matrix} \right. & \left( {{Formula}\mspace{14mu} 2} \right) \\ \left\{ \begin{matrix} {{{PosX}(k)} = {{{PosX}(n)} + {\left( {k - n} \right) \times {VelY} \times}}} \\ {\left( {{{PosY}\left( {k - 1} \right)} - B} \right)/\left( {{{PosX}\left( {k - 1} \right)} - A} \right)} \\ {{{PosY}(k)} = {{{PosY}(n)} + {VelY}}} \end{matrix} \right. & \left( {{Formula}\mspace{14mu} 3} \right) \end{matrix}$

The estimated motion of the moving body P is not limited to the uniform linear motion or the circular motion described above, and can be changed to any motion to correspond to an observation object.

The image extracting step (Step S30) extracts (trims) an image that includes the estimated position L_(n) of the moving body P from the time-series images F_(n) (n=1 to N) as an extracted image Q_(n). The extracted image Q_(n) can have any size as long as it is smaller than the size of the time-series images F_(n). The number of the time-series images F_(n) that are objects of extraction is any plural number. As this number is increased, reliability of the shape of the moving body P in a template image Q_(T) described later is increased. On the other hand, as this number is increased, a processing rate of the controller 30 is lowered. Therefore, it is desirable that the number of the time-series images used for extraction is an appropriate value determined by considering these circumstances. In the present embodiment, for convenience of description, the extracted image Q_(n) is supposed to be extracted from all the acquired N images.

The template generating step (Step S40) generates the template image Q_(T) by stacking (superposing) extracted images Q_(n) in such a manner that reference points of the moving body P coincide with each other. Each reference point is the estimated position L_(n) on the extracted image Q_(n), for example. The estimated position L_(n) has been already acquired in the position estimating step (Step S20). Therefore, the controller 30 can generate the template image Q_(T) without newly performing a calculation process. The template image Q_(T) is generated by calculating average values of pixel values of the respective extracted images Q_(n). For example, in a case where the extracted image Q_(n) has a size of M_(x)×M_(y), each pixel value T(i, j) in the template image Q_(T) is calculated by following Formula 4.

$\begin{matrix} {{T\left( {i,j} \right)} = {\frac{\sum\limits_{n = 1}^{N}{T_{Qn}\left( {i,j} \right)}}{N}\left( {{i = {1 \sim M_{x}}},{j = {1 \sim M_{y}}}} \right)}} & \left( {{Formula}\mspace{14mu} 4} \right) \end{matrix}$

Here, T_(Qn)(i, j) is a pixel value at a position (i, j) in the extracted image Q_(n). That is, an intensity distribution of pixels that represent the moving body P in the template image Q_(T) (that is, a moving body P_(T)) is an average of intensity distributions of pixels that represent the moving body P in extracted images Q₁ to Q_(N) (that is, moving bodies P₁ to P_(N)).

As illustrated in FIG. 3, the shapes of the moving bodies P₁ to P_(N) in the extracted images Q₁ to Q_(N) always change by being affected by fluctuation caused by the atmosphere or the like. The template image Q_(T) represents a state where such a shape change is suppressed. That is, the moving body P_(T) in the template image Q_(T) has an “average” shape obtained from the moving bodies P₁ to P_(N) in the extracted images Q₁ to Q_(N). For example, individual shape changes that appear in the extracted images Q₁ to Q_(N) are averaged in the template image Q_(T) to appear as “blur” that has a low intensity in the profile of the moving body P_(T). Meanwhile, the shape common to the moving bodies P₁ to P_(N) remain in the template image Q_(T) as it is. In other words, the moving body P_(T) includes a bright portion that represents the shape common to the moving bodies P₁ to P_(N) and a dark portion that represents the shape that is not common to the moving bodies P₁ to P_(N). Therefore, the moving body P_(T) can have a shape close to an original shape of the moving body P, so that the template image Q_(T) has high reliability as a template.

As described above, each pixel value in the template image Q_(T) is an average value of pixel values in the extracted images Q₁ to Q_(N). Therefore, random noise that appears in each of the extracted images Q₁ to Q_(N) is reduced in the template image Q_(T), so that an S/N ratio (a signal-to-noise ratio) is improved. FIG. 4A illustrates an intensity distribution in an extracted image Q_(n), and FIG. 4B illustrates an intensity distribution in the template image Q_(T). High-intensity portions that are observed at the centers of the intensity distributions in FIGS. 4A and 4B represent a moving body P_(n) and the moving body P_(T), respectively. As illustrated in FIG. 4A, random noise appears around the moving body P_(n). Meanwhile, as illustrated in FIG. 4B, the random noise observed in FIG. 4A is significantly reduced around the moving body P_(T). Further, because of reduction of the random noise, a variation in intensities that represent the profile of the moving body P_(T) is reduced as compared with that in intensities that represent the profile of the moving body P_(n). This reduction also contributes to high reliability as a template.

The position identifying step (Step S50) identifies the position of the moving body P in the time-series image to which the template matching is applied by performing template matching using the templates image Q_(T) for at least one of the time-series images F₁ to F_(N). This template matching is performed around the estimated position L_(n) (in a predetermined region including the estimated position L_(n)) obtained in the position estimating step (Step S20). The image to which the template matching is applied may be the extracted image Q_(n) or an image including the estimated position L_(n) which is newly extracted (trimmed) from the time-series images F_(n). In both cases, pixel values that represent the moving body P (P_(n)) and its surrounding area are the same, and therefore the same result is obtained. In the following description, it is assumed that the time-series image F_(n) is used as the image to which the template matching is applied for convenience of descriptions.

In template matching, the degree of similarity between the template image Q_(T) and the time-series image F_(n) can be evaluated by a value of zero-mean normalized cross correlation R_(ZNCC) represented by Formula 5, for example. An average pixel value in the template image Q_(T) in Formula 5 is obtained by Formula 6, and an average pixel value in the time-series image F_(n) (or an image extracted for template matching) is obtained by Formula 7.

$\begin{matrix} {R_{ZNCC} = \frac{\sum\limits_{j = 0}^{{My} - 1}{\sum\limits_{i = 0}^{{M\; x} - 1}\left\{ {\left( {{I\left( {i,j} \right)} - \overset{\_}{I}} \right)\left( {{T\left( {i,j} \right)} - \overset{\_}{T}} \right)} \right\}}}{\sqrt{\sum\limits_{j = 0}^{{Ny} - 1}{\sum\limits_{i = 0}^{{N\; x} - 1}{\left( {{I\left( {i,j} \right)} - \overset{\_}{I}} \right)^{2} \times {\sum\limits_{j = 0}^{{My} - 1}{\sum\limits_{i = 0}^{{M\; x} - 1}\left( {{T\left( {i,j} \right)} - \overset{\_}{T}} \right)^{2}}}}}}}} & \left( {{Formula}\mspace{14mu} 5} \right) \\ {\overset{\_}{T}\; = \frac{\sum\limits_{j = 0}^{{My} - 1}{\sum\limits_{i = 0}^{{M\; x} - 1}{T\left( {i,j} \right)}}}{M_{x}M_{y}}} & \left( {{Formula}\mspace{14mu} 6} \right) \\ {\overset{\_}{I}\; = \frac{\sum\limits_{j = 0}^{{Ny} - 1}{\sum\limits_{i = 0}^{{N\; x} - 1}{I\left( {i,j} \right)}}}{N_{x}N_{y}}} & \left( {{Formula}\mspace{14mu} 7} \right) \end{matrix}$

Here, M_(x) and M_(y) are the numbers of pixels that respectively represent the width and the height of the template image Q_(T). N_(x) and N_(y) are the numbers of pixels that respectively represent the width and the height of the time-series image F_(n). (i, j) is coordinates in each image (where 0≤i≤M_(x)−1, 0≤j≤M_(y)−1). T(i, j) is a pixel value (an intensity, a brightness value) in the template image Q_(T). T⁻ (in the formulae, T and a macron (a bar) above T) is an average pixel value (an average intensity, an average brightness value) in the template image Q_(T). I(i, j) is a pixel value (an intensity, a brightness value) in the time-series image F_(n). I⁻ (in the expression, I and a macron (a bar) above I) is an average pixel value (an average intensity, an average brightness value) in the time-series image F_(n).

By this template matching, the position of the moving body P_(n) in the time-series image F_(n), which is most approximated to the template image Q_(T), is identified as a proper position L_(Sn) (n=1 to N) of the moving body P in the time-series image F_(n). In other words, the estimated position L_(n) of the moving body P_(n) in the time-series image F_(n) is corrected to an identified position L_(Sn) by the template matching. The identified position L_(Sn) of the moving body P_(n) is stored in the storage unit 32 (see FIG. 1).

In the present embodiment, the template image Q_(T) is used as a comparison reference in template matching. As described above, the template image Q_(T) is generated based on a so-called stacking method (superposition method) from extracted images Q_(n) among time-series images F_(n). Therefore, even if the shape of the moving body P_(n) that is an object of comparison always changes because of an influence of fluctuation or the like, this change is hardly reflected on the shape of the moving body P_(T) in the template image Q_(T), whereas this change is strongly reflected on the shape that is common to the moving bodies P_(n). Accordingly, the template image Q_(T) has high reliability as a template, and it is possible to accurately identify the position L_(Sn) of the moving body P_(n) in the time-series image F_(n). That is, according to the present embodiment, the moving body P can be detected accurately by a stacking method.

As described above, according to the present embodiment, the position L_(Sn) of the moving body P_(n) is identified (corrected) in each of the time-series images F_(n). Therefore, by connecting the positions of the moving body P in the respective images to each other, it is possible to recognize the locus or the behavior of the moving body P that is an observation object. That is, even if a stacking method is used, fine information with regard to the position of the moving body P within an observation time can be obtained.

Each pixel value in the template image Q_(T) may be a median or a mode of intensities (pixel values) of a corresponding pixel in extracted images Q_(n) used for generation of the template image Q_(T).

The moving body observation method according to the present embodiment may include the following process. FIG. 5 is a flowchart according to a first modification of the present embodiment. As illustrated in FIG. 5, the moving body observation method may include a binarizing step (Step S60) and a center-of-gravity calculating step (Step S70). Each of the steps is performed by the moving body observation device 10 illustrated in FIG. 1.

The binarizing step (Step S60) binarizes each pixel value in the time-series images F_(n) (n=1 to N) or the extracted images Q_(n) (n=1 to N) by using a predetermined threshold. As the threshold, a value that is a predetermined multiple of a value of background noise obtained from an observation history, for example, can be applied. The binarizing step (Step S60) obtains the binarized extracted images Q_(n). The center-of-gravity calculating step (Step S70) calculates the center of gravity of the moving body P_(n) in each time-series image F_(n) or each extracted image Q_(n), which has been binarized.

The template generating step (Step S40) in the first modification uses centers of gravity obtained in the center-of-gravity calculating step (Step S70) as reference points of the moving bodies P_(n) when the extracted images Q_(n) are stacked on each other. That is, in the template generating step (Step S40), the extracted images Q_(n) are stacked on each other to make the centers of gravity coincident with each other, thereby generating the template image Q_(T). Thereafter, the position identifying step (Step S50) is performed using each of the binarized images.

FIG. 6 is a flowchart according to a second modification of the present embodiment. As illustrated in FIG. 6, the moving body observation method may include the binarizing step (Step S60) according to the first modification and a fitting step (Step S80). As in the first modification, each of the steps according to the second modification is also performed by the moving body observation device 10 illustrated in FIG. 1.

As described above, the extracted images Q_(n) that have been binarized are obtained in the binarizing step (Step S60). The fitting step (Step S80) fits a normal function to two one-dimensional distributions obtained by integrating a pixel distribution of the moving body P_(n) in each extracted image Q_(n) that has been binarized in predetermined two directions. Accordingly, a position with the maximum brightness of the moving body P_(n) (the maximum-brightness position) is calculated.

By the binarizing step (Step S60), pixels in the extracted image Q_(n) are classified into pixels with a pixel value of 1, which represent the moving body P_(n), and pixels with a pixel value of 0, which represent the background of the moving body P_(n). Therefore, in the fitting step (Step S80), a one-dimensional distribution in the X direction obtained by integrating a two-dimensional distribution of these pixels in the Y direction is calculated, and a one-dimensional distribution in the Y direction obtained by integrating that two-dimensional distribution in the X direction is also calculated. The normal function is fitted to these two one-dimensional distributions, thereby calculating the position with the maximum brightness of the moving body P_(n) (the maximum-brightness position).

The template generating step (Step S40) in the second modification uses maximum-brightness positions obtained in the fitting step (Step S80) as reference points of the moving bodies P_(n) when the extracted images Q_(n) are stacked on each other. That is, in the template generating step (Step S40), the extracted images Q_(n) are stacked on each other to make the maximum-brightness positions coincident with each other, thereby generating the template image Q_(T). Thereafter, the position identifying step (Step S50) is performed using each of the binarized images.

The present disclosure is not limited to the embodiment described above and is defined in the descriptions of the scope of claims, and the present disclosure includes all sorts of modifications with equivalent meanings and within the scope of the descriptions in the scope of claims. 

What is claimed is:
 1. A moving body observation method comprising: acquiring time-series images taken at a predetermined time interval; obtaining an estimated position of a moving body in each of the time-series images based on an estimated motion of the moving body assumed in advance; extracting images each including the estimated position of the moving body from a plurality of images of the time-series images, as extracted images; generating a template image by stacking the extracted images in such a manner that reference points of the moving body in the extracted images coincide with each other; and identifying a position of the moving body in at least one image of the time-series images by performing template matching using the template image for the at least one image of the time-series images.
 2. The moving body observation method according to claim 1, wherein each of the reference points of the moving body in stacking of the extracted images is the estimated position of the moving body.
 3. The moving body observation method according to claim 2, wherein each pixel value in the template image is set to any one of an average pixel value, a median, and a mode of a corresponding pixel in the extracted images.
 4. The moving body observation method according to claim 1, further comprising: binarizing each pixel value in the time-series images or the extracted images; and calculating a center of gravity of the moving body in each of the binarized time-series images or each of the binarized extracted images, wherein each of the reference points of the moving body in stacking of the extracted images is the center of gravity of the moving body.
 5. The moving body observation method according to claim 1, further comprising: binarizing each pixel value in the time-series images or the extracted images; and calculating a maximum-brightness position of the moving body by fitting a normal function to two one-dimensional distributions obtained by integrating a pixel distribution of the moving body in each of the binarized extracted images in predetermined two directions, wherein each of the reference points of the moving body in stacking of the extracted images is the maximum-brightness position of the moving body.
 6. The moving body observation method according to claim 1, wherein the moving body is debris on an earth orbit or a microorganism in a liquid. 